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Abstract 

I have proposed a non-Abelian and stochastic self-organized criticality model in which each 
avalanche contains one stochastic site and all remaining sites in the avalanche are deterministic 
with a constant threshold E\. Studies of avalanche structures, waves and autocorrelations, size 
moments and probability distribution functions of avalanche size, for the thresholds 4 < Ej. < 256, 
were performed. The shell-like avalanche structures, correlated waves within avalanches, complex 
size moments and probability distribution functions show multifractal scaling like the Abelian and 
deterministic BTW model despite the fact that the model is non-Abelian and stochastic with 
unbalanced relaxation rules at each stochastic site. 

PACS numbers: 45.70.Ht, 05.65. +b, 05.70.Jk, 64.60.Ak 
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I. INTRODUCTION 



Bak, Tang, and Wiesenfeld (BTW) [1] introduced a concept of self-organized criticality 
(SOC) to study dynamical systems with temporal and spatial degrees of freedom. They 
designed a simple cellular automaton with conservative and deterministic relaxation rules 
to demonstrate the SOC phenomenon. Manna (M) 2] proposed another conservative SOC 
model in which stochastic relaxation rules instead of deterministic rules were defined. 

A stable configuration (see below) in the BTW model does not depend on the order of 
relaxations so the model is Abelian On the other hand, a stable configuration in the 
M model depends on the order of relaxations and the model is thus non-Abelian. Dhar 4| 
theoretically proved the Abelian property of the M model for the case when probabilities of 
many final stable configurations are considered. 

Based on the real-space renormalization group calculations, Pietronero at al. |5| claimed 
that both deterministic [1] and stochastic [2] models belong to the same universality class, 
i.e. a small modification in the relaxation rules of the models cannot change the universality 
class. It was assumed that both models show a finite size scaling (FSS) [61. With FSS the 

n 

avalanche size, area, lifetime, and perimeter follow power laws with cutoffs [6]: 



P(x) = x~ T *F(x/L D *), (1) 

where P(x) is the probability density function of x, F is the cutoff function, and r x and 
D x are the scaling exponents. The set of scaling exponents (t x , D x ) defines the universality 
class tel. 

n n 

Avalanche structure studies [7| and numerical simulations [8[ showed that the BTW and 
M models do not belong to the same universality class. Later, Tebaldi at al. [9] found that 
avalanche size distributions in the BTW model follow multifractal scaling. They concluded 
that the avalanche size exponent r s (Eq. [T|) does not apply to the BTW model. An avalanche 
wave decomposition approach 10[ was applie d to demonstrate the different wave features 
of the BTW and M models. Karmakar at al. [12| found that models with balanced relaxation 



rules are similar to the deterministic BTW model and mode 
rules are similar to stochastic M-like models. These results 



s with unbalanced relaxation 

m 



ii 



121 ] support the idea 



that the BTW and Manna models are prototypes of different universality classes. 

The BTW and M models are well understood 0, E 9 ] . On the hand, we know very little 
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about the transition from multifractal to fractal scaling [l^l, Q]- I modified the BTW model 
[l| to study this transition. I allowed stochastic relaxations in one site in each avalanche. 
The stochastic site is located in the place where the avalanche is initiated. Avalanches are 
randomly initiated in various lattice sites so the position of stochastic site is changed at every 
new avalanche and the stochastic M-like site can visit all lattice sites for sufficiently large 
set of avalanches. An avalanche dynamics study is useful to understand the transition 
rom multifractal to fractal scaling and to test the hypothesis of precise relaxation balance 



12j . The results suggest that the model cannot belong to either BTW or M universality 
classes. 

In Sec. [TT] I introduce a non-Abelian and stochastic model. In Sec. Ill II I apply numerical 
simulations and statistical methods to find avalanche structures, autocorrelation functions, 
Hurst exponents, avalanche size moments and probability density functions of avalanche 
sizes. Section [TV] is devoted to a discussion which is followed by conclusions in Sec. [V] 

II. STOCHASTIC SELF-ORGANIZED CRITIC ALITY MODEL 

The stochastic SOC model is defined on a two dimensional (2D) lattice L x L where each 
site i has assigned two dynamical variables E(i) and E C (V) [13f| . The variable E(i) represents 
for example, energy and variable E c (i) represents a threshold at site i. All thresholds E c (i) 
are equal to the same value E\ in the interval 4 < E T C < 256. Relaxation rules are undirected, 
conservative and deterministic for all sites i with the thresholds EAi) — El. At each site 

n 

i the relaxation rules are precisely balanced [12J. Thus the sites behave as the BTW sites 



1] for the threshold E T C =4. In a stable configuration (a stationary state), all sites i follow 
E(i) < E c (i). Let us assume that from a stable configuration we iteratively select i at random 
and increase E(i) — > E(i) + 1. If an unstable configuration is reached, i.e. E(i) > E^., then 
the site i is labeled as 1m- The initial threshold at site Im^c^m) = E^ is changed to the new 
value E c (iyi) = E]} = 2 and stochastic relaxation rules |2| are assigned to this site. In each 
avalanche only one site 1m undergoes stochastic relaxations while all remaining sites relax 
as undirected, deterministic and conservative sites. All unstable sites i (including im) where 
E(i) > E c (i) undergo relaxations to reach the stable configuration E(i) < E c (i). If a stable 
configuration is reached, then the threshold E c (iu) at the site iM is set to E c (iy[) = E^ and 



deterministic relaxation rules [l] are assigned to the site iM- The site iM disappears and 
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all sites i of the lattice are BTW-like sites. Stable and unstable states are repeated many 
times. Adding of energy (E(i) — > E(i) + 1) takes place randomly thus stochastic sites 1m 
could visit all lattice sites. Thus stochastic sites introduce an annealed disorder in the initial 
deterministic model 111. 



III. RESULTS 

One stochastic site in each avalanche makes the model stochastic [yj 
All sites, except one, relax energy as Abelian sites 



Il5| . however I classified the model 



as non- Abelian. To demonstrate the non- Abelian property consider two critical sites [15] 
(Fig. [TJ with thresholds E], 1 = 2 (the stochastic site) and E\ = 4 (the deterministic site). 
One can verify (Fig. [TJ that the order of relaxations of unstable sites leads to different final 
stable configurations i.e. the OSS model is non-Abelian. 

Relaxations in all deterministic sites are precisely balanced [12] . However, one site in each 
avalanche shows unbalanced relaxations. Thus we consider the OSS model for the model 
with unbalanced relaxation rules. 

Analyzing the avalanche structures can provide an important initial information j?]. [12I \\\ 
about the nature of the SOC model. Several avalanches of the OSS model, at the threshold 
E^ = 4, have been decomposed into clusters with equal numbers of relaxations (Fig. [2]). 
Avalanche structures of two types: (i) without holes Fig. |2]Ja) and (ii) with rare holes in 
clusters Fig. Mjo) were observed. These structures are similar to shell-like structures one 
finds in the BTW model. 

The study of mathematical SOC models allows avalanches to be decomposed avalanches 



into waves J10|. This approach enables the investigation of correlation of waves within the 
avalanche jll|. I demonstrated that the OSS model is non-Abelian, thus the avalanche waves 
could be defined as in the case of the M model 11]. On the other hand, waves initiated 
by relaxations of the stochastic site (im) propagate though the lattice of BTW sites. Thus 
waves would be defined as well as waves of the BTW model. An avalanche of size s is 
decomposed into m waves with sizes Sk, where s = S^Sfe. A time-sequence of avalanche 
waves Sk is used to determine the autocorrelation function 



n( , T \ _ ( s k+tSk)L — (Sfc)| , , 

C[t ' L) ~ (4h-(s k )l ' (2) 
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FIG. 1: (color online) Order of relaxations of the critical sites shows the non-Abelian property of 
the OSS model. The circle denotes a stochastic site and the star (*) denotes a site which relaxes 
as first. In the initial state there are two critical sites Ej, 1 = 2 (stochastic) and E ! c = 4. (a) 
The stable configuration in the case when the stochastic site (Ej. 1 = 2) relaxed as first, (b) The 
stable configuration in the case when the deterministic site (E*. = 4) relaxed as first. Arrows show 
directions of energy diffusion in the stochastic site. 



where time is t = 1, 2, . . ., and the time averages are taken over 5 x 10 6 waves for lattice sizes 
L = 128,256,512,1024, 2048 and 4096. The autocorrelations C BTW (t, L) of the BTW model 
(Fig. E[a)) and C oss (t, L) of the OSS model at threshold E{ = 4 (Fig. Efb)) approach zero 
only for times and i^ff exceeding the maximum number of waves in avalanches. This 

result is a consequence of correlated waves in avalanches ll|. I have observed (Fig. |3J) that at 
a given lattice size L, the time i^ff is approximately as large as time t^^f (t^ff = t^}^). I 
note that avalanche waves in the M model are uncorrelated due to autocorrelation functions 
C M (t, L) =0fort> 1 fl|l2|. 

The autocorrelations C BTW (t, L) (Eq. [2$ of the BTW model were approximated by a 
power law /(£) ~ t~ Tc and cutoff function g(t/L Dc ) [111 ]: 
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FIG. 2: (color online) Two types of avalanche structures on a 2D lattice of size 128 x 128: (a) 
without holes in the clusters and (b) with holes in clusters. Lattice sites with the same numbers 
of relaxations are shown by the same color (rainbow pseudo-color coding). Only one site in any 
avalanche relaxes in a stochastic manner as a M site (the red area), however all remaining sites 
relax in a deterministic manner as BTW sites (E* = 4). 



C BTW (t, L) = f{t)g{t/L D <). (3) 

The autocorrelations C BTW (t, L) (Fig. [3(a)) agree well with the previous results HQ. 
However, I have found that the power law approximation f(t) ~ t~ Tc does not approximate 
the autocorrelation C BTW (t, L). I have verified (Fig. [3j that the exponential function 
f(t, L) ~ exp(— a^t) better approximates not only the autocorrelation C BTW (t, L) but also 
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FIG. 3: (color online) Log-lin plots of autocorrelation functions C BTW (t, L) and C oss (t, L) 
(threshold E{ = 4) for lattices: L = 128, 256, 512, 1024, 2048 and 4096 (from left to right). The 
autocorrelation functions C x (t, L) for given lattice sizes L and for time t > 5 are approximated by 
an exponential function C x (t, L) ~ exp(— atf) (thick lines) where decay rates are: (a) = 
0.0244 (BTW model) and (b) aggg = 0.027 (OSS model). 

C oss (t, L). 

The autocorrelations C(t, L) from Eq. [3] are used to compute the time moments 



(4) 



For lattice sizes L = 128 — 2048, threshold = 4 and for several values of q, the plots 
of log < t q > vs. logL (Fig. HI inset) were used to determine the functions crf TW (q) and 
a° (q). The plots cr BTW (q) and a° (q) exhibit linear dependence for q in the range 1.0 < 
q < 4.0. From these plots (Fig. 0} the parameters D BTW = 1.06 ±0.05, D° ss = 1.09 ±0.05, 
t btw ^_ q an( j jOSS ^_ g were determined 111 ]. 



Stochastic process are often characterized by Hurst exponents [181]. Fluctuations F(t, L) 
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FIG. 4: The plots of cr c (q) are approximated by the linear function a c (q) = piq + po- For the BTW 
model (thin line) the parameters are p\ = 1.06 ± 0.05, po = 1.06 ± 0.05 and for the OSS model 
at threshold Ej, = 4 (thick line) the parameters are p\ = 1.09 ± 0.05, and po = 1.06 ± 0.05 (this 
approximation is shown with dashed line). The inset shows log-log plots of < t q > vs. L for the 
OSS model and for the exponents q = —1, 0, 1, 2, 3 and 4 (from bottom to top). 

11]: 



F(t,L) = [(Ay(tf) L -(Ay(t))l} 1 / 2 , (5) 

are used to determine Hurst exponents where y{t) = J2k=i s k anc ^ = y{k + 1) — y(k. . 

If fluctuations F(t, L) should scale as F(t, oo) ~ t H then H is the Hurst exponent 
Two exponents H BTW = 0.89 ± 0.02 and H BTW = 1/2 (Fig. [5]) were determined for 
the BTW model and for times t < and t > . Similarly, the Hurst exponents 

H oss = Q.88±0.02 and H oss = 0.610±0.001 were determined for the OSS model (threshold 
El = 4) and for times t < i^ff and t > t^ff . Fluctuations F(t, L) were also determined 
for the other thresholds 8 < E{ < 256 (Fig. EJ. For all thresholds 4 < E[ < 256 of the OSS 
model two scaling regions of F(t, L) were identified in contrast to the fluctuation of the M 
model which exhibits the single scaling with the Hurst exponent H M = 0.53 ±0.05 (it is not 
shown in Fig. |5]). 

Moment analysis was used to study scaling properties of both BTW and OSS 

models. A property x in the FSS system obeys the scaling given by Eq. [TJ The q moments 



8 



20 



CD 



16 



14 



- s 














— BTW 




— 

J/* 
J/* 


E c =4 
-E c =8 




w/ f 


-E =16 






-E =32 






-E c =256 




'// 

/ 







4 6 
Ln(t) 



10 



FIG. 5: (color online) The fluctuations F(t, L = 4096) of the BTW and OSS model with thresholds 
4 < E\ < 256. The Hurst exponents were determined for the threshold E\ = 4: H BTW = 
0.89 ± 0.01 and H oss = 0.88 ± 0.01 , the time 1.0 < t < 50.0, H BTW = 0.503 ± 0.001 and 
H oss _ o.eio ± 0.001 for the time 1000.0 < t < 10000.0. 

of x are defined as 



•Lmax 

(x q ) = J x q P(x, L)dx ~ U^ q) , 



(6) 



where cr x (q) = (q + 1 —t x )D x and x max ~ L Dx . I calculated the moments only for avalanche 
size s. The plots log(s q ) versus logL for approximately five hundred values of the exponent 
q were used to determine the functions cr BTW (q) and crf SS (q) (Eq. [6]). These plots scale 
precisely for the BTW model (Fig. H^a)) and for the OSS model (threshold E\ = 4, Fig. 
Mh)) for all exponents 0.0 < q < 4.0 and lattice sizes 128 < L < 4096. 

The results show that af TW (1.0) = 2.03 and of ss (1.0) = 2.04 for threshold E{ = 4 
(Fig. El inset), close to the expected value crf Tiy (l) =2.0 (rj, |l6j. The function a° ss (q) 
grows faster than functions crf TW (q) and cr^(q) (Fig. [7J inset) when the exponent q > 1.0 
increases. The function daf ss (q)/dq increases with increasing q > 1.0. At q — 2.07 it 
reaches the maximum Df (2.07) = 3.17 ± 0.01 and for q > 2.07 is almost constant or 
slowly decreases. For q = 4.0, the capacity dimensions are D^f (4) = 2.76, Df TW (4) = 2.88 
Q and Df ss (4) = 3.11 ±0.02. The functions daf TW (q)/dq and da° ss (q)/dq continuously 
increase when the exponent q increases in the range 1.0 < q < 2.07. This fact is considered 
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Ln(L) 



FIG. 6: Avalanche size moments at lattice size L versus lattice size L display scaling behaviour 
for wide range of exponents q. The results are shown for (a) the BTW model (128 < L < 2048) 
and (b) the OSS model with threshold E\ = 4 (128 < L < 4096) whereas selected exponents q are 
q = 1.0, 2.0, 3.0, 3.1 and 4.0 (from bottom to top). 



to be a signature of multifractal scaling 



The capacity dimension of the Manna model 



Df (q) is constant for exponents q_> 1.0: Df(q) = 2.76 ±0.01 



12 
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12 



16|. 



161 ] , which is a typical 



property of the FSS models 

Karmakar et al. [12] claimed that if local avalanche dynamics meets criterion of a precise 
relaxation balance, then the model must show the same behaviors as the BTW model. If we 
increase the threshold E 1 > 4 and we modify the relaxation rules to meet the criterion of 
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FIG. 7: Plots of da s (q)/dq versus q for the M (dashed curves, 128 < L < 2048), BTW (thin 
line, 128 < L < 2048) and OSS models for the thresholds E\ = 4 (thick solid, 128 < L < 4096). 
The inset shows plots a s (q) versus q for the M (dashed line, 128 < L < 2048), BTW (thin line, 
128 < L < 2048) and OSS (thick solid, 128 < L < 4096) models. 



the precise relaxation balance, then the model with all deterministic sites [19| has the same 
properties as the BTW model. In this model, for thresholds 8 < E c < 256, I introduced one 
stochastic site in each avalanche to compare the behaviors of the modified SOC model (see 
Sec. [UJ) with the behaviors of the BTW model. 

Autocorrelations C(t, L) Eq. [2] were determined for the OSS model with the thresholds 
4 < E c < 256. Avalanches waves within avalanches are correlated (Fig. [H]) for all thresholds 
because for the time t < t max autocorrelations C(t, L) are greater than (C(t, L) > 0) and 
for the time t > t max the autocorrelations approach the value C(t, L) = 0.0. 

Functions da° ss (q)/dq of the OSS model for thresholds 4 < E c < 256 and function 
daf TW (q)/dq of the BTW model are shown in Fig. |9j The results show that daf (q)/dq > 
daf TW (q)/dq for thresholds 4 < E c < 256 and for exponents 1.0 < q < 4.0. A difference 
between functions da° ss (q)/dq — daf TW (q)/dq at given exponent q is observed to be higher 



than an expected experimental error 



12| 



Probability density functions P{s) of avalanche size s were determined for different lattice 
sizes L = 128, 512 and 4096 (Fig. [TO] (a)) and thresholds E* = 4, 8 and 256 (Fig. 01(b)) to 
know the impact of the lattice size L or thresholds E^ on the probability density functions 
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FIG. 8: (color online) Autocorrelations C(t, L = 4096) of the OSS model for thresholds 4 < E[ < 
256. 




FIG. 9: (color online) The plots of da s (q)/dq of the BTW and OSS (4 < E c < 256) models 
(128 < L < 4096). 

P(s). The probability density functions P(s) of the OSS model show small increases of their 
slopes for avalanches of size s < 10. 
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FIG. 10: (color online) Probability distribution functions P(s) of avalanche size s for (a) the 
constant threshold E c = 4 and increasing lattice size L = 128, 512 and 4096 and for (b) the 
constant lattice size L = 4096 and increasing thresholds E c = 4, 8 and 256. 

IV. DISCUSSION 



The avalanche structures of the OSS model (Fig. [2]) are similar to the shell-like structures 
of the BTW model[3, However, detailed analysis of these structures shows that the 
structures do not have the same properties as the structures of the BTW model. For 
example, we can see holes inside OSS structures (Fig. [2] (b)) which are not possible in the 



BTW model 



121 ] . These holes support our classification of the OSS model as a model with 
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unbalanced relaxation rules 



12|. 



Autocorrelations C(t, L) of the OSS model for thresholds 4 < E l c < 256 (Figs. Mb) 
and [H]) are C(t, L) > for times t max . This is a consequence of correlated avalanche size 
waves within avalanche Autocorrelations of the BTW model are the same as in 



the paper 



11]. However, I cannot confirm that these autocorrelations are approximated 



by the power law approximation in Eq. [31 I have observed that an exponential function 
C(t, L) ~ exp(— at) where a is a decay rate, is a better approximation of the autocorrelations 
C(t, L) then the power law approximation Eq. [31 This finding is supported by the fact that 
time-moment analysis (Fig. H|) for the BTW and OSS models leads to the time exponents 
t btw _j_ q anc j T oss _j_ q fj^g reason f or this discrepancy with the previous results llj is 
not clear and additional study is necessary. On the other hand, the maximum number of 
waves scales as t 1 ^™ ~ L an d this result confirms the previous observation [llj]. Similarly, 
for the OSS model t° s J ~ L (Fig. SJ). 

The fluctuations F(t, L) and corresponding Hurst exponents H± TW and if.f TW/ of the 
BTW model (Fig. [5]) agree well with the previous results llj,[l2|]. Fluctuations F(t, L) of the 
OSS model for thresholds 4 < E\ < 256 have two scaling regions (Fig. [5]) and corresponding 
Hurst exponents H° ss and if^ 155 . An existence of two scaling regions confirms the correlated 
avalanche size waves. 

The moments of avalanche size Fig. [7] for the BTW and M models agree well with 
previous results [III, 112]. The plots da s (q)/dq of the OSS model in Figs. [7] and [9l for 
thresholds 4 < E\ < 256 and exponents q > 1.0, are not constant as in the case of the M 
model. The increase of da s (q)/dq for q > 1.0 when the exponent q increases is considered 
for a signature of multifractal scaling [9|, Based on the moment analysis (Figs. [7] and [9]) 
for the thresholds 4 < < 256 and previous conclusions 0, [ll|, [l^] , I claim that avalanche 
size distributions show multifractal scaling. 

Probability density functions of avalanche size P{s) of the OSS model show a moderate 
increase of small avalanches s < 10 for thresholds 4 < E], < 256 (Fig. [TPl . I assume that 
these changes of probability density functions -P(s) do not influence scaling of avalanche size 
moments (Fig. [6]). 

Holes in some avalanches (Fig. ETb)) are characteristic for the models with unbalanced 

on 

relaxation rules which exhibit uncorrelated avalanche waves [3, [12| . On the other hand, the 
existence of holes in the OSS model is not sign of uncorrelated waves beacuse avalanche 
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size waves are correlated (see correlations C oss (t,L) (Fig. [3] (b) and EJ) and fluctuations 
F oss (t,L) (Fig. EJ)) exhibit two scaling regions. 

Shell-like avalanche structures j?! I12I], avalanche wave correlations C oss (t,L) 
avalanche wave fluctuations F oss (t, L) [ll| and avalanche size moments c^ 55 ^) 
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2] 
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of the OSS model support the conclusion that the OSS model shows multifractal scaling 
for thresholds 4 < E* c < 256. I can reproduce the plots daf TW (q)/dq and da° ss (q)/dq 
(Fig. [7j) for thresholds 4 < E* < 256. Comparison of the functions daf TW (q)/dq with the 
previous results of the BTW model and undirected model [3] shows that these functions 
collapse to a single function. I demonstrated that the OSS model shows multifractal scaling 
but the functions daf ss (q)/dq, for thresholds 4 < E T C < 256, are not identical with the 
function daf TW (q)/dq of the BTW model (Figs. [TJ The differences between functions 
daf TW (q)/dq and da® (q)/dq at given exponent q are larger than the method error 12]. 
Based on these evidences, I conclude that OSS and BTW models belong to the multifractal 
universality class. However, the models do not have identical properties despite the fact 
that in the OSS model only one site in each avalanche undergoes stochastic relaxations as 
the M site. 



V. CONCLUSION 

The OSS model has been developed to study properties of the inhomogeneous sand pile 
model [l3] at very low densities of M sites 



20|. Based on traditional classification schemes 



[12| , one can expect that the model will belong to the M universality class. However, I have 
demonstrated that one stochastic M site in each avalanche is not enough to change multifrac- 
tal scaling of the model to the FSS (Fig. [7]). The OSS model is stochastic, non-Abelian with 
unbalanced relaxation rules [ijj], the classification schema implies that the model belongs 
to the M universality class. However, the OSS models exhibits correlated avalanche waves 
and multifractal scaling which is not allowed for the models in the M universality class. The 
OSS model exhibits multifractal scaling and the classification scheme implies that the model 
belongs to the BTW universality class, however the model is stochastic, non-Abelian, and 
has an unbalanced relaxation rule [lj| , thus it cannot belong to the BTW class. I think that 
it could be more convenient to consider multifractal or FSS scaling as a main criterion for 
model classification to solve this paradox. Then the models which show multifractal scaling 



15 



(Fig. [7J BTW and OSS models) despite the fact that they are not identical (see Sec. 
could belong to the multifractal universality class. Models that show FSS scaling Eq. [I] 
(M model) could belong to the FSS universality class. I have analyzed another SOC model 
where the results support this classification scheme [19I ] . 
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This work was supported by the Slovak Research and Development Agency under contract 
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